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Fronts in randomly advected and heterogeneous media and 
nonuniversality of Burgers turbulence: Theory and numerics 
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A recently established mathematical equivalence — between weakly perturbed Huygens fronts (e.g., 
flames in weak turbulence or geometrical-optics wave fronts in slightly nonuniform media) and the 
inviscid limit of white-noise-driven Burgers turbulence — motivates theoretical and numerical esti- 
mates of Burgers-turbulence properties for specific types of white-in-time forcing. Existing math- 
ematical relations between Burgers turbulence and the statistical mechanics of directed polymers, 
allowing use of the replica method, are exploited to obtain systematic upper bounds on the Burgers 
energy density, corresponding to the ground-state binding energy of the directed polymer and the 
speedup of the Huygens front. The results are complementary to previous studies of both Burgers 
turbulence and directed polymers, which have focused on universal scaling properties instead of 
forcing-dependent parameters. The upper-bound formula can be heuristically understood in terms 
of renormalization of a different kind from that previously used in combustion models, and also shows 
that the burning velocity of an idealized turbulent flame does not diverge with increasing Reynolds 
number at fixed turbulence intensity, a conclusion that applies even to strong turbulence. Numeri- 
cal simulations of the one-dimensional inviscid Burgers equation using a Lagrangian finite-element 
method confirm that the theoretical upper bounds are sharp within about 15% for various forcing 
spectra (corresponding to various two-dimensional random media). These computations provide a 
new quantitative test of the replica method. The inferred nonuniversality (spectrum dependence) 
of the front speedup is of direct importance for combustion modeling. 

PACS numbers: 05.20.Jj, 36.20.Ey, 42. 15. Dp, 47.70.Pq 



I. INTRODUCTION 

In the geometrical optics of a "quenched" medium with 
static spatial variations in refractive index, or in the com- 
bustion of a heterogeneous solid propellant, each spatial 
point is characterized by a local speed v(x) at which a 
wave front or flame sheet can advance normal to itself in 
accordance with Huygens' principle A key problem 
is to determine the time needed for an influence (light 
or combustion) to propagate from one region of space 
to another; Fermat's principle asserts that the propaga- 
tion effectively occurs along all possible trajectories that 
obey the local speed v(x), with the first-arriving trajec- 
tory giving the desired answer. 

A preliminary connection with statistical mechanics is 
suggested by writing the travel time for a spatial path C 



t(C)= f ds ^ = (v- 1 )l(C)+ f dsa(x), (1) 

Jc n x ) Jc 



is a spatial aver- 
-i\ 



where 1(C) is the path length, (v ) 
age over the medium, and a — v^ 1 — {v~ M ) is a zero- 
mean field. We can interpret C as a material curve — a 
polymer string — and t(C) as its potential energy, with 
(v^ 1 ) being the constant tension and a an external po- 
tential. The polymer is infinitely flexible, since there is 
no curvature energy in Eq. {1}. And the polymer is "di- 
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rected" because each endpoint is confined to a different 
region of space. Fermat's principle requires us to find the 
minimum-energy configuration (classical ground state) of 
this directed polymer. 

So far this is just a complicated optimization prob- 
lem, but it can be made more explicitly statistical in two 
ways. First, in many cases of interest, a is a homoge- 
neous random field with specified statistics, and we are 
satisfied with calculating the ensemble-averaged propa- 
gation time or minimum energy (which, for a very long 
polymer, nearly equals the minimum energy in a single 
realization of a). Second, even within a given realization 
of the potential, it may be convenient to consider the 
canonical ensemble of directed-polymer configurations at 
finite temperature, calculate the thermodynamic free en- 
ergy, and then take the zero-temperature limit to obtain 
the ground-state energy. The most powerful theoretical 
methods arise from combining both types of ensemble. 

This paper focuses on the case in which a is a homoge- 
neous and isotropic (but not necessarily Gaussian) ran- 
dom field that is small in magnitude compared to (v^ 1 ). 
In this weak-fluctuation limit, there is in fact a more gen- 
eral and useful relation between propagation and directed 
polymers, one that applies even to random media with 
time dependence and advection, such as turbulent fluids. 
In this paper, we focus on the quenched-medium inter- 
pretation of the analysis, but relate our approach to tur- 
bulent combustion and other applications when useful. 
As was recently shown [2j, the propagation of Huygens 
fronts in a broad class of weakly random d-dimensional 
media reduces to the (d — l)-dimensional inviscid Burg- 
ers equation with white-in-time forcing. Furthermore, 
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the statistical properties of the inviscid Burgers dynamics 
are equivalent to those of viscous "Burgers turbulence" 
(with the same forcing) in the limit of vanishing viscosity 
v Q. The white-noise-driven viscous Burgers equation, 
in turn, has a known relation to the canonical ensemble 
for a near-straight directed polymer (at a temperature 
proportional to v) in a ti-dimensional static random po- 
tential that is white in one direction [J] . This is precisely 
the kind of random potential commonly assumed in the- 
oretical directed-polymer studies [1, IE 0] • 

In Sec. |TT]we systematically discuss some known rela- 
tions among different theoretical representations of our 
problem, and indicate how the recent results can be 
linked with this framework to allow calculation of the 
speedup of Huygens fronts. In Sec. IIIII we describe three 
rigorous "monotonicity properties" obeyed by random 
Huygens propagation, which provide not only checks on 
our later results but also further implications from them. 
In Sec. lIVI we present the replica method in a form suited 
to our problem, and derive explicit upper bounds on the 
front speedup. In Sec. [V] we review existing numerical 
simulations that provide relevant speedup values for com- 
parison, and describe our new higher-precision simula- 
tion results. A summary and discussion are presented in 
Sec. EH 



II. RELATIONS AMONG THEORIES 

A. KPZ, Burgers, and Huygens 

The Kardar-Parisi-Zhang (KPZ) equation for interface 
growth [8| is well known as a unifying model for diverse 
phenomena. The equation can be written 



dh 



v\7 2 , h 



(2) 



where /i(i,xjj is a fluctuation in the height of an inter- 
face as a function of d— 1 transverse dimensions and time, 
v > is a surface tension that smooths the interface, and 
r)(t, xj_) is a zero- mean external perturbation (not nec- 
essarily the white noise assumed by KPZ). Equation ([2|) 
is equivalent to the (d — l)-dimensional forced viscous 
Burgers equation 



9w 
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Vjjw = vV\vf - V±r] 



(3) 



for the velocity field w = — Vj_/i, and also to the (d — 
l)-dimensional imaginary-time Schrodinger equation for 
the "wave function" exp(h/2v), with potential —rj and 
Planck's constant 2v. The solution of this Schrodinger 
equation is given by the Feynman path integral Q 



exp 



2v 



Mo,y(o)) 



2v 



Vy{u) exp 
y(t)=x_L v -'' 

t du[i|y'(«)| 3 -/ / (</.y(«)); ). (4) 
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The simplest connection between these ideas and Huy- 
gens propagation is to assume a quenched medium with 
weak refractive-index fluctuations given by a smooth 
function cr(xu , Xj_), and to write the position of a nearly 
flat front as the unperturbed position plus a small dis- 
placement: 



(5) 



in units where (v~ x ) = 1. We can approximate the in- 
cremental speed T](t, xj_) in Eq. ([2]) as the negative of 
the index fluctuation at the unperturbed location, i.e., 
T](t, xj_) ~ — <r(t, xjj. Equation @ then governs the 
dynamics of h, with the nonlinear term describing the 
leading-order effect of tilted propagation, and the Lapla- 
cian term smoothing the cusps (or Burgers shocks) that 
would otherwise develop. (We work to leading order 
in the fluctuation a and the tilt Vj_/i, neglecting their 
higher-order and joint effects.) The parameter v, which 
corresponds to the Markstein length in premixed flamelet 
combustion (loj . contributes a stabilizing dependence of 
front speed on local curvature, with concave regions prop- 
agating faster than convex regions. Pure Huygens prop- 
agation is obtained as a "viscosity solution" in the v — ► 
limit [Hill. 

A useful alternative description is based on a field 
To(x||,xj_) giving the time at which the front reaches a 
given point. From Eq. ([5]), we see that to leading order 



T (x||,xj_) = X|| - ft,(x||,xj_). 



(6) 



Thus, in the weak-fluctuation framework, we are free to 
write X|| in place of t as an independent variable and 
interpret — h as the time deviation. With this perspec- 
tive, the path integral ^ is, upon suitable normaliza- 
tion, the canonical partition function for the simple "Fer- 
mat's principle" directed polymer described in the Intro- 
duction. With weak fluctuations, only paths that are 
nearly aligned with the xm -direction are competitive in 
travel time (energy). Such a path, described by a func- 
tion y(u) so that x^ = y(x||), has arc-length element 
ds = du [1 + ||y'(w)| 2 ] and travel time 



i(y(u)) = / ds[l + a(u >y {u))} 



du[l + hy'(u)\ 2 +a(uMu))}. (7) 



Then, assuming a flat initial front with h — 0, Eq. (|U) 
asserts that 
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d U [I|y'( U )| 2 + ( 7( U ,y( U )); 



(8) 



With the travel time as the energy, the right-hand side 
is the partition function at temperature 2j/ for a poly- 
mer with one end constrained to the initial front and the 
other end fixed at (xu,x±). Consequently, 7o(x||,x_l) is 
the free energy of this polymer. As expected, Eq. (jHJ) in- 
dicates that To approaches the absolute minimum travel 
time (ground-state energy) in the Huygens-propagation 
limit v — > 0. 

If a is a homogeneous random field, we expect the free 
energy To to scale linearly with the longitudinal extent 
X\\ of the polymer in the thermodynamic limit xm — > oo. 
Since the ground-state energy in the absence of fluctua- 
tions would be just x\\, we define the binding energy per 
unit length as 



A = lim 1 - 



T (X||,X_l) 



lim 



>0, 



(9) 

which is independent of by homogeneity, and is pos- 
itive because the nonlinear term in Eq. @ makes h in- 
crease on average. This nonlinear term equals \w 1 , and 
so A can also be described as the steady-state energy 
density of the Burgers fluid. [The other two terms on the 
right-hand side of Eq. ([2|) average to zero because h is 
statistically homogeneous in xj_ and r\ is a centered per- 
turbation.] The effect of the weak fluctuations cr, whose 
rms value we denote by e <C 1, is to renormalize the over- 
all front speed upward to = lim a: .._ >0O x\\/Tq = 1 + A. 
Our analysis applies to the asymptotic limit e — * and is 
expected to be accurate when the dimensionless param- 
eter e is small; previous numerical work [Isj] discussed 
in Sec. IV Al indicates that the asymptotic scaling holds 
within a few percent for e < 0.1. 



B. White-noise reduction 

For media of a given structure, with fluctuations re- 
lated by overall rescaling, it is clear that A — > as e — > 0; 
but the form of this dependence is subtle. If v is fixed, 
then for sufficiently small e, Eq. ^ can be linearized; 
the lowest-order solution for h is proportional to e and 
thus the nonlinear term (producing the secular growth 
of h measured by A) is proportional to e 2 . This is the 
weak-perturbation scenario normally associated with the 
KPZ equation, corresponding to a laminar solution of 
the viscous Burgers equation (J3J). On the other hand, 
if we take v — > at finite e to obtain Huygens propa- 
gation, and only then take e — > 0, the linearization is 



invalid (because the Burgers flow is fully turbulent) and 
the scaling of A with e is not immediately obvious. Sub- 
ject to mild conditions on the medium structure, we have 
shown mathematically 0] that A cx e 4 / 3 in this regime, 



confirming a previous conjecture [14|. Here we provide a 
somewhat more physical argument based on dimensional 
analysis. 

Let us denote by oeu and a± some measures of the lon- 
gitudinal and lateral correlation lengths of a, and con- 
sider the behavior of the inviscid Burgers equation upon 
varying these quantities and e, while keeping fixed all 
other details of the statistics of a. We can perform di- 
mensional analysis on the KPZ equation (JSJ, with rj = —a 
and v = 0, by assigning the conventional dimensions from 
the Burgers-fluid interpretation: 

[to] = [V_l/i] = LT-\ [h] = L 2 T"\ [ay] = T, 

[a x ] = L, M = [e] = L 2 T- 2 , [A] = L 2 T~ 2 . (10) 

(These dimensions are unusual from the viewpoint of 
propagation because the longitudinal and lateral direc- 
tions are treated differently.) The most general dimen- 
sionless combination of input parameters is a function of 



q — a 2 a ± 2 e. Accordingly, 



A 



(11) 



where we have taken advantage of the freedom to choose 
any combination of parameters with the dimensions of A 
to multiply Q(q). (We could equally well have chosen e 



or a,, a 



The reason for our choice is that if we take a 

-1/2 







while scaling e cx a,, , so that a becomes white noise 
in the longitudinal (time) direction, then the coefficient 
of Q(q) is unaffected. Assuming that the inviscid Burg- 
ers equation is well behaved with white-in-time forcing, 
so that A remains finite, we can partially constrain the 
function Q. The fact that e — > oo as an — > is of no con- 
cern here, because the sole formally dimensionless mea- 
sure of the strength of e is q, which is approaching zero. 
We conclude that Q(q) has a finite limit as q — > 0. Now, 
keeping the correlation lengths fixed but taking e — > 
(which also gives q — > 0), we see that A cx e 4 / 3 . 

Our mathematical analysis [2j demonstrated this re- 
sult more seamlessly, starting from the exact equations 
of Huygens propagation (rather than from the KPZ equa- 
tion), extracting a factor e 4 / 3 by a rescaling of vari- 
ables, and directly obtaining the white-noise-driven in- 
viscid Burgers equation as e — > 0. (To make the argu- 
ment used in this paper more complete, in Appendix [X] 
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we provide a mathematical justification for starting from 
the KPZ equation.) Furthermore, we showed that the 
effect of an advecting velocity field u of order e — along 
with possible time dependence of a and u on a natu- 
ral timescale of order e _1 — is given as e — ► by simply 
adding —um to a and considering the medium as quenched 
at the initial time. That is, both the component and 
the time dependence of the medium become irrelevant. 

If we reintroduce the viscosity (Markstein length) v, 
with conventional dimensions L T , there is a new di- 
mensionless parameter that remains fixed in the white- 
noise limit: the Burgers-fluid Reynolds number Res ~ 
(aiio^e 2 ) 1 / 3 ^ -1 . To maintain Huygens propagation (cor- 
responding to very large Ree) as e — > 0, we must decrease 
v at least as fast as e 2 / 3 . In our previous analysis, we 
cited mathematical results [H, [l5[ establishing that the 
white-noise-driven steady state of the inviscid Burgers 
equation exists, and is approached by that of the viscous 
Burgers equation as Ree — > oo. Thus there is no subtlety 
with interchange of the limits e — * and Ree —> oo, as 
there was with e — > and v — > in Eq. @ . Previously [2| 
we denoted by v what we here call ve~ 2 ^ oc Re^ 1 . Let 
us emphasize that Ree is distinct from the hydrodynamic 
Reynolds number ReNS f° r weak advection by a turbulent 

I 



Here we have observed that the —ik\\X\\ term is negligi- 
ble because the correlation function with the argument 
£-2/3,2, vanishes as e — > for all nonzero xu , and we have 
then rescaled the integration variable xn. The spectrum 
is white in the longitudinal direction because there is no 
dependence on ku. The e~ 2 factor simply normalizes the 
original a correlation function to unit variance. 

Thus the desired prefactor of e 4 / 3 in the Huygens-front 
speedup is the value of A obtained using the path integral 
(TJI in the v — > limit, with the spectrum \Y1\ for r\. In 
other words, the prefactor is the ground-state binding 
energy per unit length for the "directed polymer" with 
energy 

H = J du[\\y>{u)\ 2 - V {u,y{u))]. (13) 

(The ground-state energy would be zero in the absence 
of fluctuations, and so the binding energy is simply the 
negative of the free energy.) Although this is the system 
most commonly referred to as a directed polymer in a 
random potential, we observe that the white noise rj is 
in no sense small, and thus 5|y'(w)| 2 is not an accurate 
approximation to an incremental arc length. The energy 



Navier-Stokes flow. In such a flow, we have an ~ a± ~ L 

and ReNs ~ £ ^ns (where we have now returned to units 
in which (v^ 1 ) = 1, so that L — T and e is dimensionless). 
Assuming near-equality between the flame parameter v 
(arising from thermal diffusivity) and the hydrodynamic 
viscosity v^s, as is valid for gaseous combustion, we have 
Re B ~ e 2/3 Lt/ N g > Re NS for e < 1. Thus the Burg- 
ers turbulence is more fully developed than the hydrody- 
namic turbulence, and the Huygens limit Rcb — > oo is 
not incompatible with the physics of gaseous flames. 

The quantitative implications of the white-noise reduc- 
tion follow from Eq. (fTTj) . We consider altering a given 
physical field a (rms value £ < 1) by multiplying the 
correlation length a\\ by e 2 ^ 3 (compressing in the longitu- 
dinal direction) and multiplying a everywhere by e~ 4 / 3 . 
The new field then has an rms value e = e -1 / 3 . The pa- 
rameter q for the new field is unchanged from the original, 
and so A is multiplied by e -4 / 3 , based on the coefficient 
of Q(q). Thus we have removed the e 4 / 3 dependence, 
and the new A is the prefactor. In the e — > limit, the 
spectrum of the noise 77(^11, xx) = — £~ 4 / 3 <j(e~ 2 / 3 :E|| , xj_), 
obtained by altering the original field a as described, be- 
comes 



(12) 



(j 1 3[) is to be considered on its own terms, separate from 
the original intuition about a flexible string with ten- 
sion. In particular, to use a supposedly more accurate 
square-root expression for arc length in Eq. (|13p would 
be incorrect for our purposes. Remarkably, the ideal- 
izations conventionally adopted to render the physical 
directed polymer more tractable (quadratic expansion of 
arc length, white- noise potential) are precisely justified in 
the problem of weakly perturbed Huygens propagation. 



III. MONOTONICITY PROPERTIES 

A. Dependence on perturbation spectrum 

It is intuitively reasonable that a more vigorous forcing 
of the Burgers equation should result in a higher steady- 
state energy density, and correspondingly that a random 
medium with greater fluctuations should result in faster 
propagation. More precisely, let us consider the v — * 
limit of the path integral (jHJ) for two different homoge- 
neous weakly random media, labeled 1 and 2, and assume 
that at every wave vector the spectrum of o\ is no greater 



L>(fc||,k ± ) = J da; || d rf_1 xx exp(-ifc||a;|| - ik_L • xj_) e~ 8/3 (er(0, 0) cr(e" 2/3 X|| , x x )) 
= J dx\\ d d_1 x_L exp(-ik^ ■ xjj e -2 (cr(0, 0) er(a;||,xj_)}. 
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than that of 02 ■ Because the spectra of independent ran- 
dom fields are additive, we can obtain a field £ with the 
same spectrum as 02 by defining 

Z = a 1 +p, (14) 

where p is a zero-mean Gaussian random field (indepen- 
dent of <7i) with an everywhere nonnegative spectrum 
given by subtraction. According to Eq. (fT2|) . as e — * 0, 
the speedup is governed entirely by the spectrum (or 
equivalently the two-point correlation function) of the 
medium fluctuations, and so £ gives the same speedup 

as (72. 

Now we observe that the fastest path from the initial 
front to a given point in the field 01, which we denote 
yi(w) (with travel time t\), can be considered as a possi- 
ble path in the field £. There, its travel time is different 
by an amount 

r x \\ 

1=1 dup(u,y!(u)). (15) 
Jo 

The fastest path in the field £ must have a travel time 
no greater than t\ + I. The ensemble average over £ can 
then be performed in two steps, thanks to the statistical 
independence of p and <j\. First, upon averaging over 
p for a given realization of a\, the mean of I is zero 
and so the "partial average" first-arrival time for £ is 
no greater than t±. Second, upon averaging over a±, the 
mean of t\ is (obviously) the average first-arrival time for 
cr±. Consequently, the "complete average" first-arrival 
time for £ (or (72) is no greater, and the speedup no 
smaller, than for a±. 

This monotonicity property does not immediately re- 
late the prefactors A of e 4 / 3 for different media, because 
two different spectra with the same total power e 2 can- 
not obey the assumed inequality. But relations can be 
obtained by considering media with e differing as little 
as possible such that one spectrum is bounded by the 
other, and then correcting the result with the known e 4 / 3 
scaling. The necessary discrepancy in e (which weakens 
the prefactor relation) can be reduced by performing an 
optimally chosen spatial rescaling of one of the media, 
which does not affect its Huygens-front speedup. 

Table Q] lists the correlation functions and spectra of 
four homogeneous isotropic random media in d = 2 di- 
mensions, which we will use as examples throughout this 
paper. Each medium is parametrized by the rms fluctua- 
tion e and a length scale a. We note that the white-noise 
spectrum (fT"2|) is just the d-dimensional Fourier transform 
of the normalized correlation function at the wave vec- 
tor k = (0,k^). For isotropic media, this is a function 
D(|kj_|) of exactly the same form as the Fourier trans- 
form D(k) at a general wave vector. The unnormalized 
spectrum to be used with the monotonicity property is 
e 2 D(k). Table HT1 gives the results of numerically optimiz- 
ing the relations between pairs of spectra to constrain the 
speedup prefactors A. The omitted pairs are not useful 
because one spectrum dominates as k — ► and the other 
as k — ^ 00. 



B. Dependence on spatial dimension 

A second monotonicity property applies when a given 
medium can be viewed as a "slice" through a higher- 
dimensional medium M, i.e., as the restriction of the field 
a to a (hyper)plane perpendicular to the initial front in 
M. All paths in the slice are also paths in M with the 
same travel time, and the longitudinal coordinate xii is 
the same in both media. It follows that the first-arrival 
time in M is no greater, and the speedup no smaller, 
than in the slice. 

For homogeneous isotropic random media, the two- 
point correlation function of a in the slice is exactly the 
same function of distance as in M. Because of the dif- 
ference in dimensionality, however, its Fourier transform 
(i.e., the spectrum) is generally different. Of course, a 
correlation function is realizable only if the associated 
spectrum is nonnegative. We see that if a given function 
of distance is a realizable correlation function in M, it 
is automatically realizable in lower dimensions, and thus 
the spectrum (although different) remains nonnegative. 
But in higher dimensions, the spectrum may develop neg- 
ative values and the correlation function may no longer 
be realizable. 

The most powerful application of this monotonicity 
property is to a correlation function that is realizable 
in all dimensions, such as the Gaussian, whose Fourier 
transform is always another Gaussian function. In this 
case, the speedup prefactor must be a nondecreasing 
function of d for all d > 2. 



C. Dependence on laminar flame speed 

The special case of Huygens propagation with advec- 
tion by a random velocity field, but with a fixed front- 
advancement speed in the local comoving frame, is widely 
studied as a model of turbulent premixed combustion 
0, There, the fixed local speed is called the lami- 
nar flame speed ul, and the overall statistically steady 
propagation rate is called the turbulent burning velocity 
ut- An important question in combustion modeling is 
the dependence of ut on ul and on the statistics of the 
advecting flow field (such as rms velocity vl), assumed 
for simplicity to be given in advance rather than dynam- 
ically affected by the flame. We have shown [2[ that 
for u 1 Jul <C 1 (weak turbulence), this idealized com- 
bustion problem is equivalent to weakly random Huy- 
gens propagation in a particular quenched medium de- 
termined by the flow statistics (specifically, the two-point 
spatial correlation function). Thus the methods of this 
paper, phrased for convenience in terms of quenched me- 
dia, apply also to flames in weak turbulence. 

Moreover, a simple monotonicity property allows lim- 
ited conclusions even about the opposite limit m'/ul S> 1 
(strong turbulence, the more important case in practice) . 
Namely, if all flow statistics are held fixed, ut must be 
a nondecreasing function of ml. This holds realization 
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TABLE I: Correlation functions and normalized spectra of example two-dimensional random media (all homogeneous and 
isotropic) . 



Medium 


<a(0)a(r)> 




D(k) 


Gaussian (G) 


e 2 exp(— r 2 /a 2 ) 




2 { 1 2 7 2\ 

7ra exp( — |a k ) 


Modified Gaussian (MG) 


e 2 (f-r 2 /a 2 )exp(- 


-7a 2 ) 


±7ra 4 fc 2 exp(-± a 2 fc 2 ) 


Exponential (E) 


e 2 exp(— r/a) 




27ra 2 (f + a 2 k 2 )-^ 2 


Modified exponential (ME) 


e 2 (f-ir 2 /a 2 )exp( 


-r/a) 


Tva 4 k 2 (7 + 2a 2 k 2 )(l + a 2 k 2 )- 7 ^ 2 



TABLE II: Upper bounds (UB) on relative speedup from spec- 
tral monotonicity. 



Medium 1 


Medium 2 


£i/e2 


a 1 /a 2 


UB on Ai/A 2 


G 


E 


0.723 


1.414 


1.541 


MG 


E 


0.596 


2.000 


1.993 


MG 


ME 


0.659 


1.944 


1.743 


ME 


E 


0.904 


1.029 


1.144 



by realization: If we consider two fronts (initially coin- 
cident) propagating independently in the same flow, the 
front with smaller ul can never get ahead of the other 
front anywhere, because at the location and time of any 
such overtaking, it would have to be advancing faster 
relative to the flow than its rival. 

The connection between weak and strong turbulence 
is then as follows: For given flow statistics, the turbu- 
lent burning velocity ut is no greater for very small ul 
(a strong-turbulence problem) than for very large ul (a 
weak-turbulence problem). Numerically this relation is 
useless, since for very small ul we expect ut ~ u', with a 
coefficient that remains unknown. We can only say that 
ut < ul, where ul S> u' is a laminar flame speed large 
enough that the corresponding turbulent burning veloc- 
ity is essentially ul- But the relation is useful in ruling 
out the possibility that ut diverges to infinity at fixed 
v! as some other flow parameter is varied. For this to 
happen in strong turbulence [v! /ul S> 1), it would also 
have to happen in weak turbulence (u! /ul <C 1), i.e., the 
speedup prefactor A would have to diverge. This can be 
ruled out in a given case by a suitable upper bound on 
A. 



IV. REPLICA ANALYSIS 
A. Overview of the replica method 

Originally developed in the study of discrete statistical- 
mechanical systems such as spin glasses [lg, [l7|, [l8| , the 
replica method is a powerful tool for theoretical investi- 
gation of the thermodynamics of disordered systems. It 
exploits the interplay of two kinds of randomness: the ex- 
ternal stochastic parameters that determine the energy 
"landscape" in which a system finds itself, and the finite- 



temperature thermal fluctuations of the system configu- 
ration within that landscape. Although we are here inter- 
ested mainly in ground-state properties of the stochastic 
landscape, it is beneficial to consider a finite temperature 
and later take it to zero. 

To describe the replica method in the context of our 
problem [1, lj, HH, let us denote by Z{rj) the partition 
function for the directed polymer described by Eq. (fl~3|) 
with a particular realization of the white-noise potential 
?7(x|| , xj_). We are working with a polymer of a given lon- 
gitudinal extent at a finite temperature r. The most gen- 
eral description of the disordered thermodynamics would 
consist of the probability density function (pdf) of Z{rj) 
resulting from the given statistics of 77. This pdf would 
be difficult to obtain directly, but its moments can be 
simplified considerably. 

For any positive integer n, we have 

Z( v ) n = J V{y c (u)} e X p(-i^#(y a (u))\ (16) 

^ a— 1 ' 

a product of independent path integrals over n "replicas" 
of the polymer; here T>{y c (u)} = n"=i ^yc( u )- Because 
H is linear in 77, and because r\ has Gaussian statistics, 
we can average Eq. (|16j) using the identity 

(expC) =exp(i(C 2 )), (17) 

valid for any zero-mean Gaussian variable £. We identify 
the Gaussian variable appearing in the exponential in Eq. 
CU), 

1 ™ r 

( = -22 du v(u,y a (u)), (18) 

T a=l 

and compute 

(C 2 ) = - E / dudu'(r,(u,y a (u))r,(u',y b (u'))) 




(19) 

We have used the relation (rj(u, y) t](u' , y')) = 6(u — 
u') V(\y - y'|), where 

V(y) ^ J |^ exp(ik J ..y)U(|k J .|). (20) 
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Upon combining Eq. (fTTf with the part of H that is 
independent of 77, we obtain 

W) = / Z>{yc(«)} fflq>(-iiT n ({y c (u)})), (21) 

which has the form of a single partition function for n 
replicas with "energy" 

J ^ a=l 

1 - \ 

-57 E ^(ly-(«)-yt(«)l))- ( 22 ) 



.6=1 



This H n , unlike H, contains no random potential, but the 
price is that we now have several interacting polymers (a 
different number of them for each different moment of Z 
we wish to calculate). 

Much as in Sec. Ill A) the new partition function (Z n ) 
can be identified as a Feynman path integral for the 
n-particle imaginary-time Schrodinger equation, with 
Planck's constant t and with the potential — V/2t acting 
on every pair of particles (including self-interaction). The 
potential function is now static and deterministic, unlike 
the time-dependent random potential —77 in the original 
Schrodinger equation. The wave function ^(yi, . . . , y n ) 
is governed by the Hamiltonian operator 
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1 

0,5=1 



(23) 



The thermodynamic limit of an infinitely long polymer 
has a simple interpretation: By evolving the imaginary- 
time Schrodinger equation 



W 1 . 
a — = n-nV 

axu t 



(24) 



for an infinite time, we project any initial wave function 
onto the quantum ground state. The path integral for 
(Z n ) over a sufficiently long time Z11 is dominated by a 
term proportional to exp[— E g (n) xu/t], where E g (n) is 
the ground-state energy of the n-particle system, because 
the contributions of the other (higher) energy eigenvalues 
are asymptotically negligible. 

What we actually wish to calculate is not a moment 
(Z n ) but the averaged free energy — r(lnZ). Unfortu- 
nately, In Z cannot be expanded in a Taylor series about 
Z = 0, and to make progress we must introduce a pecu- 
liar feature of the replica method. The needed average is 
expressed by the identity 



(InZ) = lim -> '- . 

n— >0 n 



(25) 



and it is assumed that (Z n ), computed as above for pos- 
itive integer n, can be analytically continued to n near 



zero. Then the asymptotic (large Xu) behavior of the 
averaged free energy is 



ry, y ^ \~ E Q (") X \\ I A ~ 1 

-T{mZ) = —t hm 

n^O n 

— x\\ lim , 



(26) 



n— >0 n 

and so the polymer's binding energy per unit length is 



A = - lim 



E g (n) 



n—>0 n 



(27) 



B. Variational treatment 

We desire a method for estimating E g (n) that is 
demonstrably valid for any positive integer n and that 
also (unlike, e.g., numerical solution of the n-particle 
Schrodinger equation) can be formally generalized to 
noninteger n. A very useful choice is the variational 
method, which is based on the observation that the ex- 
pectation value (ip\H. n \'ip) of the n-particle Hamiltonian 
in an arbitrary quantum state |^>) is an upper bound on 
E g (n). To obtain as tight a bound as possible, this ex- 
pectation value is minimized over a convenient family of 
"trial" wave functions, in the hope that some of them are 
close to the true ground state. 

The variational method is valid for positive integer 
n, where 7i n is a Hermitian operator on a well-defined 
Hilbert space. The wave function associates a number 
with each configuration of n particles in d — 1 dimen- 
sions, but because the Hamiltonian is translation invari- 
ant, the center of mass separates and the nontrivial part 
of the wave function depends on n — 1 vectors. Further- 
more, the wave function is subject to one normalization 
constraint, (ip\ip) — 1. Thus the number of degrees of 
freedom in the wave function can be written 



(d_i)(„_i) 



1. 



(28) 



When n = 1, for example, we have / = (no 
degrees of freedom), as expected because there is a 
unique translation-invariant one-particle wave function, 
the eigenstate of zero momentum. This variational esti- 
mate is automatically the exact ground state of the free- 
particle Hamiltonian Tli. 

For n — > 0, we have / = — 1, and the replica method 
relies on the following nonrigorous argument: If E g (n) 
is the minimum of (if)\T-C n \ip) with respect to a separate 
variation in each degree of freedom of \ip) around the 
true ground state, then when the degrees of freedom are 
themselves negative in number, any conceivable variation 
of \ip) will result in a decrease of (^p\H n \ip}. Consequently, 
we maximize (ip\7i n \ip) among trial wave functions and 
thereby obtain a lower bound on E g (n). The results of 
applying this strategy to spin glasses have been verified 
by rigorous methods [20l [2lj . and there is no evidence 
that the corresponding results for directed polymers are 
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invalid. In fact, the numerical simulations in Sec. [V] pro- 
vide a successful quantitative test of this application of 
the replica method. 

In Appendix [B] we analyze a commonly used family of 
trial wave functions, parametrized by a function z(u), for 
which (t/j\H n \i/j) can be expressed analytically in n. The 
result is 



= \r 2 {d-l)n 



du 



; A(u) 



2r 



where 



A(u) = huz(u) — | / dv z(v) 



B (z) 



(27r) d - 



B o (0)+ / rf Uj B (z(u)) , (29) 



(30) 



^exp(-i|k x | 2 z)D(|k x |). (31) 



As n — > 0, a meaningful wave function requires z(u) to 
be nonnegative and nonincreasing for < u < 1. If we 
define 



r(r, z) = — lim 



1-2 



r^(d-l) 



o u 2 A(u) 



+ ^ [Bo(0) - J duB (z{u))j , (32) 

then the polymer's binding energy per unit length at tem- 
perature r is bounded above by T(r, z) for any nonneg- 
ative, nonincreasing function z(u). Thus the Huygens 
prefactor A obeys 



A < lim T(t,z t ), 



where we anticipate based on previous results Q that a 
useful, finite bound will require a r-dependent choice of 
z(u). 



C. Explicit replica bounds 



The replica variational treatment has been applied in 
some detail to the directed-polymer problem, focusing 
on the case of the Gaussian medium [7]. In that work, 
all values of temperature were considered, and the goal 
was a general qualitative understanding of the polymer's 
behavior, rather than a calculation of its binding energy. 
Here, by concentrating on the binding energy in the zero- 
temperature limit, we are able to extend the previous 
results to obtain explicit bounds on A, not only for the 
Gaussian medium but for arbitrary spectra. ( "Gaussian" 
refers to a type of spectrum, defined for d — 2 in Table 
HI As discussed in Sec. Ill a quite separate feature of 
Huygens propagation is the equivalence of weak isotropic 
perturbations, not necessarily Gaussian in pdf, to the 
directed polymer's white-noise perturbations, which are 
automatically Gaussian in pdf and described completely 
by their spectrum.) 

We start with a derivation of the equation for station- 
arity of T(t, z) and its consequences, along the lines of 
the previous work 0. From Eq. (fB6|) . we compute 



SA(v) 
5z(u) 



huS(u — v) + | 9{u — v), 



(34) 



(33) where 8(x) equals 1 for x > and for x < 0, and thus 



Sz(u) 



= -±r 2 (d-l) 
= -±r 2 (d-l) 



dv 



o v 2 A(v) 
1 

uA(m) 2 ' J v 2 A(v) 



2 lU S(u -v) + 8{u ~v)]-^ B' (z(u)) 



du 



+ — Si («(«)). 



(35) 



We use the notation 

qd-U 



= J ^pi ( 2 |k X | 2 ) p exp(-I|k x | 2 z)C(|k x |); 

(36) 

note that each B p (z) is a positive, decreasing function 
for z > 0. 



In the absence of constraints on z(u), the variational 
optimum would satisfy ST/Sz — 0. Inserting Eq. (|35p 
and differentiating with respect to u, we obtain 



1-2 



T\d-1) 



A'(u) 



uA(u) 3 2t 
or, using Eq. (|B13j) . 



B 2 (z(u))z'{u) = Q, (37) 



z'(u) 



l (d-l) B 2 (z(u)) 



16A(u) 3 



•It 



= 0. 



(38) 
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Regions where z'(u) — (for which infinitesimal vari- 
ations could violate the nonincreasing constraint) still 
obey Eq. (|38|) . Provided z(u) > everywhere (so that 
the nonnegative constraint has no local effect), Eq. ([55]) 
is a necessary condition for an optimum. 

The type of solution obtained for the d — 2 Gaussian 
medium p} has the second factor in Eq. |38|) equal to 
zero for < u < u c , while z(u) = z(u c ) for u c < u < 1 
so that the first factor z'(u) equals zero there. Thus we 
have 



A(u) = \t ( 



A(u) = \ z{u c ) 



d-l 
B 2 (z(u)) 



1/3 



(0 < u < u c ), (39) 
(u c < u < 1). (40) 



From the assumed continuity of z(u) at u c , the continuity 
of A(u) follows, and so 



\ z(u c ) = A(u c ) = 



As t — > 0, this gives 



l ( d- 1 



\B 2 (z(u c )) 



1/3 



(41) 



z(u c ) = T 



d-l 



1/3 



Also, by differentiating Eq. (f39|) . we find 



(42) 



A'(u) 



1 (d-l) 1 / 3 
6 B 2 (z(u))^ 



B 3 (z(u)) z'(u) (0<u<it c ), 

(43) 



or, assuming z'{u) ^ and using Eq. (1B13|) . 

i (d-l) 1 / 3 



" 3T B 2 (z(u)) 4 /3 



B 3 (z(u)) {0<u<u c ). (44) 



This constitutes an implicit solution for z(u). We would 
like to make use of this solution as far as possible for 
arbitrary spectra, without adopting particular forms of 



Bp. 



Let us define 



H{z) 



B 3 (z) 

B 2 { z y/3- 



(45) 



The simplest case occurs when /i(z) is a decreasing func- 
tion for all z > and lim^oo u(z) = 0. Then a single- 
valued, decreasing function z(u) for < u < u c is defined 
by 

u = \r(d- 1) 1/3 m0) [z(u c ) < z < oo], (46) 

where z(u c ) is given by Eq. (|4"2")) . We now attempt to sub- 
stitute this trial solution into Eq. |3^|) . Upon integration 



by parts, the first term of T(t, z) becomes 

1 f 1 du f Uc du 



r x = |^(d-i) 



A(uc) J Uc " 2 Jo u2 Hu) 



\r 2 {d-l) 



l/Uc - 1 

A(w c ) 



A'(u) 1 
«A(m) 2 uA(u) 



kr 2 (d-l) 



A(u c ) 



du . , /„ + hm 



A(u) 2 u^omA(w) 



(47) 



To evaluate the limit, we observe that u — » corresponds 
to z — > oo. If £>(fc) ~ k q as A: — > 0, then for the total 
power 



d d k 



(27T) 



(48) 



to be finite, as required to define the parameter e of the 
random medium, we must have q + d > 0. The resulting 
behavior of B p , from Eq. (|3"6"|) , is 

B p (z) ~ z -P-(9+ d -!)/2 (z -►«>). (49) 

Using Eqs. ([Ml) and flU, we find 



wA(tt) oc 



B 2 {zf/z 



,(?+<*) /3 



oo (z — > oo), (50) 



and so the limit in Eq. (|47p is zero. Furthermore, from 
Eq. (|1T|), the term oc t 2 /A(u c ) in Eq. (j47]) scales with r 
and vanishes as r — ► 0. Under a change of variable to z, 
the remaining integral gives 



(51) 



r 1 = \(d-i) 1 ' a / dzB 2 (zf/\ 

The second term of T(r, z) becomes 

T 2 = 1-(bo{0)-B (z(u c )) J du 



It 



du Bq(z(u)) 



B o (0) - {1 - u c ) B (z(u c )) 
+ f dzB (z)u'(z)\ (52) 

where u(z) is given by Eq. ()46jl . Integration by parts 
then yields 



r 2 = i fs (0) -S (z(«c))+ / dzBj 

2t V Aw 



(z)u(z) 



(53) 
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because lim*_»oo B$(z) u(z) — 0. Recognizing a difference 
quotient, which as r — > becomes a derivative oc Bq(0), 
we obtain 



r 2 = is^o) 



d-l 
B^fi) 

+ Hd-i) 1/3 



1/3 



2 („ c ) B 2 (zY' A 



In both Eqs. (|5T|) and (|54|) . the lower limits of the inte- 
grals can be taken to zero by Eq. (|4"2")1 , since there is no 
remaining singular dependence on r. Also, because 

d B l( z ) _ o r .N2/3 , 1 gl(gWg) 



Eq. (|54j) simplifies upon a further integration by parts to 

r 2 = i(d-i) 1 / 3 / dzB 2 ( z ) 2 / 3 

Jo 

+ \{d~lf' 3 l im _^£) . (56) 



The limit vanishes because 



B2( ^. 2 -(^)/3^ ( ^oc). (57) 
Thus we obtain the replica bound 

POO 

A<r 1 +r 2 = |(d-l) 1 / 3 / dzB 2 {z) 2 ' 3 (58) 

Jo 

on the prefactor of the Huygens-front speedup. 

A remarkable renormalization interpretation of Eq. 
is seen by rewriting it as 

A < f (d - l) 1 / 3 [ X - [z 3 ' 2 B 2 {z)f/\ (59) 
Jo z 



Note that 

z 3 ' 2 B 2 {z) 



d d ~ l \s. 



(2tt; 



^i|kJ 4 z 3 / 2 ex P H|kJ 2 z),D(|kJ) 



11/2 r(|d) 

2 r(i(d-i)) 



> 2 1 
d d k 

(27r) d 



(fc 2 z) 3/2 exp(-±fc 2 z) D(fc), 



(60) 



where we have inserted factors to compensate increasing 
the integration from d — 1 to d dimensions. Because the 
integrand becomes small for k -C z~ x l 2 or k ^§> z -1 / 2 , 
Eq. (|60[) represents the power contained in a finite wave- 
number band around k ~ z -1 / 2 of width 5k ~ z -1 / 2 . In 
Eq. (fB"9")) . this spectral band power (analogous to e 2 ) is 
raised to the | power (consistent with e 4 / 3 scaling) and 



then integrated over all logarithmic length scales (dz/z). 
This suggests a stepwise process in which, starting from 
the smallest length scales, each order-unity spectral band 
has the same qualitative effect as if acting alone: It renor- 
malizes the front propagation speed, with the new ef- 
fective speed (turbulent burning velocity in combustion) 
providing the raw input (laminar flame speed) for the 
next larger-length-scale band. In the weak-perturbation 
limit of Huygens propagation, because all these renor- 
malization contributions are very small, they combine 
additively as displayed in Eq. (|59"j) . These conclusions 
are compared to existing concepts of front-speed renor- 
malization in Sec. I VII 



The assumptions about fi(z) stated below Eq. ([45)) hold 
for the two-dimensional Gaussian and exponential media 
in Table |U As a result, bounds on their speedup pref ac- 
tors can be obtained from Eq. d5T 



An < 



38/3^1/3 



16 

A E < 2.038. 



1.714, 



(61) 
(62) 



For the exponential medium, due to the long spectral 
tail at high wave number, B p (0) is divergent for p > 
1. Although such quantities appeared in the derivation, 
the exponential medium can be approached by a limiting 
process to make the expressions well-defined. The end 
result, Eq. (|58|) . is finite and was evaluated directly by 
numerical integration to obtain Eq. (|62[) . 

We now consider violations of the previous assump- 
tions about fJ,(z). If fJ,(z) is decreasing for all z > but 
lim z _ >nn jJ>{z) — /i* > 0, then we define z(u) by Eq. (|4^)) 



for ±r(d - l) 1 / 3 ^ 



< u < u c , and define z(u) = 00 



for < u < u*. It follows that A(u) = 00 for < u < u*, 
and so both integrals in Eq. (|32p receive nonzero contri- 
butions only from u* < u < 1, since i3o(oo) = 0. The 
change of variable to z and the integrations by parts pro- 
ceed as before, and the result is unchanged. 

On the other hand, if (J,(z) is not an everywhere de- 
creasing function, let [0, z*] be the largest interval from 
zero on which it is decreasing (possibly z* = 0). We 
define = \r(d — l) 1 / 3 ^i(z st ) and again take z{u) — 
A(it) = 00 for < u < «*. Then the z integrals extend 
only up to z* , and several boundary terms from integra- 
tion by parts no longer vanish. Specifically, we find 



Ti = 



1^2 



(d-l) 



1 



M* A(u*) 



+ i(d-l)V3/ rfzi? 2 (z) 2 / 3 , 
Jo 



1 



r 2 = — B (z,)«. + |(d-l) 



1/3 



2r 



ffi(z«) 

B 2 ( z *y/ 3 



+ |( d _ 1)V3/ dzB2{z) 2/3 i 

Jo 



(63) 



(64) 



11 



and thus 

A < Ti +r 2 

, d „i/a ( 3 M**?'* x B 
~ [d ij B s (z*) +6 i 

i 



(z») S 3 (z*) 



S 2 (^) 1 / 3 



B 2 (z*) 4 / 3 
dzB 2 (z) 2 / 3 ). (65) 



This bound applies to the two remaining media in Table 
UJ For the modified Gaussian medium, fx(z) is in fact an 
increasing function for z > 0, so we take z* = 0, giving 



< 



3597T 1 / 3 



. , , , -l 1.599. (66) 
28/334/351/37 v > 

For the modified exponential medium, numerical evalua- 
tion shows that fj,(z) is decreasing only for < z < z* ~ 
7.492, giving 



A ME < 1-943. 



(67) 



The reason /i(z) ultimately increases for these media is 
that D(Jfe) ~ k 2 as k -> 0. From Eq. flUD with q = 2 and 
d = 2, we find that /i(z) ~ z 1 / 6 as 2 — > 00. By contrast, 
two-dimensional media with D(k) ~ fc° as A; — ► have 
~ z -1 / 6 as 2: — * 00. 
It is possible to obtain a slightly better bound on Amg- 
Because z* = for this medium, no real use has been 
made of the solution (|44p . The trial functions we have 
constructed are simply piecewise constant, of the form 



z(u) 
z(u) 



00, 

Zn, 



A(u) 
A( U ) 



00 
1 „ 



(0 < u < u c ), (68) 
(u c < u < 1). (69) 



We can discard the earlier motivation and consider arbi- 
trary permissible values of u c and z c . Upon variational 
optimization, Eqs. ([68|) and |69|) are known as the one- 
step solution because only one level of replica symmetry 
breaking is needed (K — 1, mi = u c ). This solution was 
previously discussed for the directed polymer with arbi- 
trary perturbation spectrum Q and will now be derived 
in our notation. Equation (|32[) becomes 



r(r,^) = ir 2 (d-l) 



1 



2 



1 



+ ^-[B Q (0)-(l-u c )B o (z c )]. (70) 
It 



Stationarity with respect to u c and z c gives 

1 , B (z c 



u*z c 



2r 



0, (71) 



-\r 2 {d-l) 
As t — * 0, the solution is 

'd-l 



1 



l + (l- Uc )^M = 0. (72) 

Z c ZT 



2 

d-l 



1/3 



Bx(0) 



1/3 



Bo(0) 2 / 3 ' 

1/3 go(O) 1 / 3 
^(O) 2 / 3 ' 



(73) 
(74) 



and substituting into Eq. (|70|) shows that 



. / d- 1 
2 \ 2 



1/3 



J B (0) 1/3 -B 1 (0) 1/3 . (75) 



For the two-dimensional modified Gaussian medium, 
Eq. (|75"|) gives a tighter bound than Eq. ([6^)l . 



34/3^1/3 
A MG < ; 1-585, 



(76) 



as expected because we have optimized over a new fam- 
ily that includes our previous trial solution. For the 
Gaussian medium, we obtain a looser bound than be- 
fore (A G < 3tt 1 / 3 /2 4 / 3 ~ 1.744); for the exponential 
and modified exponential media, since -Bi(O) diverges, 
the one-step upper bound is infinite and uninformative. 



D. Implications of the replica results 

We now discuss the replica results in light of the 
monotonicity properties of Sec. IIIIl The finite-band- 
renormalization interpretation described below Eq. (|58[) 
indicates that media with broader spectra (on a loga- 
rithmic wave-number scale) should have larger bounds 
on A. This is because the spectral power is more widely 
dispersed among bands, giving a lower amount per band 
before each is raised to the | power, and x 2 ^ 3 decreases 
more slowly than becomes small. [The limit 

x — > corresponds to a spectrum D{k) cx k~ d with power 
spread equally over many orders of magnitude in wave 
number.] Indeed, the replica bounds on A are larger for 
the "multiscale" media E and ME, which have long spec- 
tral tails at high wave number, than for the "single-scale" 
media G and MG, whose spectra fall off very rapidly 
at high wave number. Furthermore, each "modified" 
medium, exhibiting a suppression of low wave numbers, 
has a smaller replica bound than the medium from which 
it was derived. Consequently, the spectral-monotonicity 
bounds given in Table HT1 though valid, are not particu- 
larly sharp. Those bounds involve scaling down the am- 
plitude of the narrower spectrum (medium 1) until it fits 
under the broader spectrum (medium 2). The resulting 
upper bound on A1/A2 is necessarily greater than unity, 
whereas the true value is expected to be less than unity if 
these inferences based on the replica bounds are accurate, 
a hypothesis that is tested numerically in Sec. [V] 

The dependence on spatial dimension discussed in 
Sec. lIIIBl is confirmed by the replica bounds for Gaussian 
media for various d. Although in Sec. IIV CI we numeri- 
cally computed the replica bounds only for certain two- 
dimensional media, we emphasize that Eqs. (|58p. (|65[) . 
and (f75)) are valid for all d > 2 under the stated as- 
sumptions and definitions. Gaussian media are the most 
straightforward to consider for arbitrary d because the 
spectrum and the spatial correlation function can both 
retain a Gaussian form. (The other media in Table U 
are inherently two-dimensional; extending them to d > 2 
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would be a matter of definition and would require qual- 
itatively changing either the correlation function or the 
spectrum, the former possibly departing from monotonic- 
ity, the latter jeopardizing realizability. ) The normalized 
Gaussian spectrum is the d-dimensional Fourier trans- 
form of exp(— r 2 /a 2 ), i.e., 

D(k) = J d d r cxp(— ik ■ r) exp ^— —^j 

= ^/Vexp(-±a 2 /fc 2 ). (77) 

For d > 3, we find that (J,(z) is a nondecreasing function 
(much as for the modified Gaussian medium in d = 2), 
and we are driven to the one-step solution. Equation (f75|) 
gives 

A<|^(d-l) 2 / 3 ^ 1.744(d-l) 2 / 3 , (78) 

consistent with A being an increasing function of d. [The 
extra factor (d — l) 1 / 3 comes from the d-dependence of 
Bp.] The one-step replica solution for Gaussian media 
with d > 3 was previously discussed from the perspective 
of white- noise-driven Burgers turbulence Q • 

Finally, the replica bounds, in combination with the 
link between weak and strong random advcction in 
Sec. IIII C[ have an important implication for idealized 
turbulent combustion. Weak-turbulence bounds follow 
from the relation between weakly random quenched and 
advected media mentioned in Sec. Ill Bl Specifically, an 
isotropic incompressible random flow, with kinetic energy 
per unit wave number E(k), is equivalent in speedup to 
an isotropic quenched medium with unnormalized spec- 
trum 

This relation leads to e 2 = it' 2 /(d — 1) in terms of the 
mean square velocity it' 2 = 2 / °° dk E(k) — as expected 
because at each relevant wave vector k (with ku = 0), 
the velocity fluctuations arc distributed over d — 1 direc- 
tions transverse to k, only one of which (the xn direction) 
contributes to the speedup. 

The Kolmogorov spectrum of Navier-Stokes turbu- 
lence in the limit of infinite Re has E(k) ~ k d+1 for 
k — ► and E(k) ~ /c~ 5 / 3 for k — > oo. We can take 
Re = u'L/i/ns — > oo by reducing z/ns or — since the 
speedup is insensitive to spatial rescaling — by increas- 
ing the integral scale L at fixed Ujsss- This Kolmogorov 
spectrum is qualitatively similar to that of the modified 
exponential medium. The effective D(k) for turbulence 
is proportional to k for k — > (matching medium ME) 
and to fc~ 2 / 3 ~ d for k — > oo [i.e., a one-dimensional spec- 
trum fc d_1 D(k) ~ fc~ 5 / 3 , versus k~ 2 for medium ME]. 
Just as with medium ME, B p (0) diverges in Re = oo tur- 
bulence for p > 1, so the one-step bound on the flame 
speedup is uninformative. But Eq. (I65|) provides a finite 
upper bound, because the Kolmogorov spectrum gives 



-62(2) ~ z 7 / 6 as z — > (for all d) and thus the integral 
of Bi{z) 1 ^ converges. The value of z* is finite because 
fx(z) decreases for small z \p(z) = —B^IB^ ~ z~ n / 18 ] 
and then increases for large z [fi(z) ~ ^t^- 1 )/ 6 from Eq. 
(|49|) ]. We conclude that the weak-turbulence speedup is 
finite even for Re = 00. Since the turbulent burning ve- 
locity ut is a nondecreasing function of the laminar flame 
speed Ul , it follows that ut remains finite for Re = 00 in 
the case of strong turbulence (ul — * for a given flow, 
i.e., fixed u'). The importance of this result for combus- 
tion modeling will be discussed in a future publication. 

V. NUMERICAL TESTS 
A. Existing simulations 

Here we review the available data that quantitatively 
describe the behavior of a relevant system (one of the 
class of equivalent problems discussed in Sec.[TTJ and can 
be compared directly with our analytical results. These 
data are from numerical simulations with d = 2: cither 
two-dimensional weakly random Huygens propagation or 
one-dimensional white-noise-driven Burgers turbulence. 
Existing experimental front-propagation results are not 
sufficiently reliable for comparison due to the difficulty 
of approaching all the required idealizations (pure Huy- 
gens propagation, weak perturbations, unbounded statis- 
tically homogeneous medium). Even d = 2 experiments 
that appear to confirm the e 4 / 3 speedup scaling [22j do 
not correspond to the simple white-noise reduction of 
Sec. Ill 131 because the medium is artificially constructed 
from statistically independent patches on a regular grid 
aligned with the propagation direction, producing long- 
range correlations in the medium structure. Numerical 
simulations of three-dimensional Huygens propagation 
[2 13 also yield e 4 / 3 scaling but involve similar long- 
range correlations in addition to transverse anisotropy. 

Comparisons with two of our example spectra can be 
made for existing simulations of Huygens propagation in 
two-dimensional isotropic quen ched media, motivated by 
applications in seismology [13j. Medium G and medium 
ME (there called simply "exponential") are synthesized 
as Gaussian random fields, with e ranging from 0.005 
to 0.1, and travel times are computed by an algorithm 
based on Huygens' principle. Plots show the expected 
transient growth of the speedup, with a significant (but 
not yet complete) leveling-off at the longest propagation 
distances. Thus the results obtained should underesti- 
mate the steady-state speedup, and definitely be below 
the replica bounds. The power laws reported from fit- 
ting the speedup are 0.0026(100e) 133 (med ium G) and 
0.0035(100e) 126 (medium ME). Adjusting the results to 
an e 4 / 3 law based on a central value e = 0.02 to obtain 
the best estimate of the prefactor, we find 

A G > 1.2, (80) 
A ME > 1.5, (81) 
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where the inequalities reflect the incomplete equilibra- 
tion. Indeed, these lower bounds are consistent with, 
and reasonably close to, the replica upper bounds ((6T|) 
and 1(57)1 ■ 

Useful results for medium MG are available from high- 
resolution numerical simulations [23|, [24| of the one- 
dimensional viscous Burgers equation with white-in-time 
forcing at Ree ~ 10 4 (very close to the inviscid limit), 
where Ree is the Burgers-fluid Reynolds number defined 
in Sec. Ill Bl The spatial forcing spectrum corresponds 
to medium MG, but, as with other Burgers simulations 
focusing on universal features like small-scale structure 
functions and velocity pdf tails, the key nonuniversal pa- 
rameters (forcing amplitude, energy density) are reported 
only roughly. Using raw simulation data [251 ] . however, 
we determine the steady-state energy density 



-,{w 2 ) = 1.60(20), 



(82) 



where the la statistical uncertainty is estimated by di- 
viding the data into three segments, and is substantial 
because only a few "large-eddy turnover times" are sim- 
ulated in a steady state. (This is an appropriate tradeoff 
for simulations focusing on small-scale features and thus 
requiring high resolution.) Equation ((52")) applies for a 
forcing spectrum [24], [25[ with a = 4x 10 4 and 



dk± 



\k\D{\k x \) = 5 x 10~ 13 = 



32 



(83) 



(exact values), whereas the normalized MG spectrum in 
Table |T] has this integral equal to 157r 1 / 2 /a 3 . Correcting 
Eq. ((82]) with the factor (15tt 1 / 2 /32) 2 / 3 , we obtain 



A MG = 1.42(17). 



(84) 



This result suggests that the replica bound ((75)) is valid, 
and if valid, it is seen to be fairly sharp. 




1.4 



MG 



G ME 



FIG. 1: Values of A for media in Table U (note reordering). 
Shaded bars: region excluded by replica bounds. Symbols: 
numerical results (error bars indicate combined statistical and 
systematic uncertainty) . 



The numerical results are plotted, along with the 
replica bounds of Sec. IIVCI in Fig. [1] Our numerical re- 
sults are consistent with, but substantially more precise 
than, the existing simulations described in Sec. IV Al The 
replica bounds are seen to be not only valid but also sharp 
within about 15%. Furthermore, the relative order of A 
among the media agrees with that of the replica bounds, 
supporting the validity of the finite-band-renormalization 
picture discussed in Sec. IIV1 The significant variation of 
A among media, and the close agreement with replica 
bounds, suggest that the replica formulas are useful and 
accurate also for the practically important case of Huy- 
gens propagation in d = 3, although no reliable data for 
comparison are known. 



VI. DISCUSSION 



B. New simulations 

To obtain high-precision values for the speedup pre- 
factor in our four example media, we have developed a 
geometric algorithm for numerical evolution of the one- 
dimensional inviscid KPZ equation 

dh , / dh\ 2 , . , . 

«=HaO +rKM) (85) 

with white-in-time forcing. This can be simulated more 
efficiently than the original propagation problem, be- 
cause for very small e, the white-noise process in t re- 
flects a longitudinal distance scale of front evolution that 
is much longer than the correlation length of the medium 
[l|. To define the problem precisely, we assume stan- 
dard periodic boundary conditions on a lateral domain 
< x < L, with bulk properties recovered in the limit 
L — > oo. Details of the numerical method are given in 
Appendix 



Motivated by the problem of Huygens-front propaga- 
tion in isotropic random media, which reduces to pre- 
viously studied white-noise systems (Burgers turbulence 
and directed polymers) in the weak-perturbation limit, 
we have performed a systematic study to obtain quanti- 
tative information about the front speedup via the prefac- 
tor A of the already established e 4 / 3 scaling [2] . The pre- 
factor A corresponds to the energy density of the Burg- 
ers fluid and the binding energy of the directed poly- 
mer, making these "toy models" directly applicable to 
the propagation problem. The latter, though also ide- 
alized, is physically more realistic because, e.g., white 
noise is not assumed. We have extended the variational 
analysis based on the replica method — previously applied 
to directed polymers [f| 0, @] and then to Burgers tur- 
bulence [H — with a specific focus on the value of A in 
the zero-temperature or inviscid limit, corresponding to 
Huygens propagation. This analysis has been found suf- 
ficiently tractable to yield explicit upper bounds on A 
for arbitrary perturbation spectra D(k), subject to the 
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previously identified conditions 2] for reduction of the 
propagation problem to white noise. 

Let us summarize how the numerical value of the 
replica bound on A can be obtained for a d-dimensional 
random medium with a particular spectrum D(k), ei- 
ther specified analytically or determined from an exper- 
iment or simulation to sufficient precision to perform 
the required computations. [A turbulent energy spec- 
trum E(k) can be converted to an equivalent quenched- 
medium spectrum using Eq. (|79p .] First we obtain the 
functions B p {z) from Eq. ([36]) for p = 0, 1, 2, 3, and de- 
fine /i(z) by Eq. ([45]) . A general replica bound formula is 
then Eq. (f65[) . where > is such that u(z) is a decreas- 
ing function for < z < z*. If fi(z) is decreasing for all 
z > 0, we can take z„ = oo and use a simpler formula, Eq. 
(|58| . An alternative bound applicable in all cases, which 
may be better or worse (or even completely uninforma- 
tive), is the "one-step" result, Eq. (|75[) . Heuristically, the 
one-step bound tends to dominate for single-scale media 
in high dimensions or in which low wave numbers are 
suppressed (the two are related because the "volume" of 
low wave numbers has a k d factor). 

The replica results are particularly interesting in light 
of rigorous properties of random Huygens propagation 
that have been deduced from general arguments. The 
dependence of the replica bounds on the form of the spec- 
trum (confirmed by numerical results) indicates that the 
rigorous bound on relative speedup derived in Sec. IIII Al 
is not usefully sharp when applied to spectra of differ- 
ent shapes. The replica bounds, at least for one class 
of media, are consistent with the fact that a lower- 
dimensional "slice" through a medium has a smaller (or 
equal) speedup due to elimination of some possible paths. 
For randomly advected Huygens propagation (considered 
as an idealization of turbulent combustion), monotonic- 
ity with respect to the laminar flame speed ul, in con- 
junction with the finite replica bound for weak Re = oo 
turbulence, precludes a divergent turbulent burning ve- 
locity ut in strong Re = oo turbulence (ul ~ * 0). 

The key qualitative insight obtained from the analytic 
form of the replica bounds is the concept of finite-band 
renormalization. The picture of a progressively coarse- 
grained medium with an upwardly renormalized propa- 
gation speed has been previously used in turbulent com- 
bustion [26|, [27| , but it was assumed that the renormal- 
ization is purely local in wave number, i.e., that the effect 
on the renormalized speed ur of eliminating an arbitrar- 
ily narrow high-wave-number band depends only on the 
spectral power <j) in that band. On dimensional grounds, 
then, the change in renormalized speed is 

( <t> \ r/2 

Su R oc ur I —2 J , (86) 

WrJ 

where we assume that the dependence on <f> is a power 
law. If the spectrum consisted only of the band in ques- 
tion, then Eq. ([56]) would have to reproduce the weak- 
perturbation speedup scaling (now known to be r = |). 



Equation (|86[) can be rewritten as a simple additive renor- 
malization, 

5(u R ) oc (ff/ 2 . (87) 

It was then argued [13] that in turbulent combustion, 
the effect of eliminating many such bands in succession, 
covering an entire spectrum, would be to increase ur 
from ul to ut, with a cumulative renormalization of u r R 
proportional to (it' 2 ) r / 2 (where u' 2 is the total spectral 
power), giving 

u r T - u r L oc u' r . (88) 

This formula has the encouraging feature that ut oc vl 
for ul — > 0, as expected. 

We observe, however, that the use of arbitrarily narrow 
wave- number bands is incompatible with r = |. If the 
spectrum is divided into M bands of equal power <f> = 
u' 2 /M, then the effect of M renormalizations by Eq. ([57| 
is 

u r T - ul cx M f-j , (89) 

which has a finite limit as M — > 00 only if r = 
2. This exponent value, corresponding to an e 2 weak- 
perturbation speedup, was in fact suggested by an ear- 
lier field-theoretic renormalization analysis based on in- 
finitesimal wave-number bands [26j , which is widely used 
as a model of ut- That analysis, besides having the 
wrong weak-turbulence scaling, predicts that ut depends 
only on ul and u' but not on the form of the spectrum, 
in contradiction to the nonuniversality (spectrum depen- 
dence) seen in our analytical and numerical results. 

The unsuitability of infinitesimal bands for analyzing 
Huygens propagation is seen not only formally but also 
physically. The rationale for stepwise renormalization 
is that the front reaches a steady state with respect to 
small-scale perturbations, thereby determining the ef- 
fective speed of a coarse-grained front that responds to 
larger-scale perturbations. This picture is literally appli- 
cable if the perturbations exist on two widely separated 
scales. For continuous spectra, such renormalization is 
approximately justified if the bands are wide enough to 
give significant scale separation, but narrow enough to 
be roughly monochromatic so that spectral shape is not 
an issue within each band. Thus it is not surprising that 
the replica bounds involve the power in an order-unity 
band, Eq. ([50]). 

Front-speed renormalization has here been placed on a 
sounder footing by means of the replica method, but only 
in the weak-perturbation limit. It is tempting to conjec- 
ture that a relation like Eq. ([88]) with r = | may still hold 

beyond that limit, with u^ 3 — u 4 ^ 3 given by a quantity 
characterizing the random advection. This would imply 
that the weak-turbulence speedup prefactor A can be 
used to determine the strong-turbulence value of ut by 
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taking ul — > 0. Such a relation, however, is incompati- 
ble with the expectation that uj< in strong turbulence de- 
pends not only on the spectrum (or two-point spatial cor- 
relation function), which completely determines A, but 
also on other flow properties including time dependence 
and higher moments. Time dependence clearly can af- 
fect ut because, e.g., if the flow correlation time goes 
to zero at fixed u' and ul, then the turbulent diffusivity 
vanishes and the effect of advection disappears. Depen- 
dence on only the two-point spatial correlation function 
is an asymptotic result of the central limit theorem for 
u'/ul — * and does not apply beyond that regime Q. 
Thus a relation like Eq. ([88)) can hold only for a restricted 
class of flows, if at all. More generally, it is not yet clear 
whether the finite-band renormalization concept is useful 
beyond the weak-perturbation limit. 



Because the accuracy of the replica bounds on A is 
not rigorously established in general, it is important to 
seek independent validation. To this end, we have pre- 
sented the results of high-precision numerical simulations 
of the one-dimensional inviscid KPZ-Burgers equation 
with white-in-time forcing, which corresponds to two- 
dimensional Huygens propagation. The numerical results 
are within about 15% of the replica bounds for four ex- 
ample media. Although the replica method has been vali- 
dated rigorously for spin glasses [2(| Hl| , the present work 
is the only quantitative test of directed-polymer replica 
bounds known to us. 



High-precision simulations in a greater number of di- 
mensions would be very costly. Well-controlled experi- 
ments would be an alternative way to validate the replica 
results for three-dimensional propagation. One specula- 
tive possibility is based on reinterpreting the function 
To appearing in the eikonal equation (|A1|) as an electro- 
static potential. If a heterogeneous ferroelectric medium 
can be constructed in which the electric field — VTo at 
each point has a frozen magnitude cx 1 + but an un- 
constrained direction, then one face of the medium can 
be grounded (To = 0), which determines the electric- 
field direction throughout the medium, and the potential 
can be measured along the opposite face to determine 
the "speedup." (A small dissipative term cx V 2 T is 
needed in the eikonal equation to regulate singularities 
in accordance with Huygens' principle, and so the local 
charge density must affect the mechanism that freezes the 
local electric-field magnitude.) By whatever technique, 
confirmation of the sharpness of the replica bounds for 
three-dimensional propagation under idealized conditions 
would establish these bounds as an appropriate starting 
point and limiting case for more complex and realistic en- 
gineering models, such as are needed in turbulent com- 
bustion. Other physical applications were noted previ- 
ously 0. 
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APPENDIX A: APPLICABILITY OF THE KPZ 
EQUATION 

In Sec. [II] we assumed that the KPZ equation ade- 
quately describes the propagation of an initially flat front 
in a quenched medium with weak random fluctuations. 
Here, the justification of this assumption is explained. An 
exact equation for Huygens propagation in a quenched 
medium is the eikonal equation [2] 

\VT \=l + a, (Al) 

where To(x) is the arrival time at a point x and cr(x) is 
the refractive-index fluctuation. If we define h{x\\ , x^) = 
£11 — To(xii,xj_) in accordance with Eq. @ and assume 
that the overall propagation is in the +a;||-direction, then 
the eikonal equation can be written 

= 1 - V(l + <7) 2 -|V^| 2 . (A2) 

In an initial interval of xm where |V_l/i| 2 is small com- 
pared to typical values of a, the right-hand side of Eq. 
(|A2j) equals —a to leading order. Thus the tilt Vi/i, 
initially zero, grows with xy at a rate proportional to 
the amplitude of a (measured by the rms fluctuation 
e <C 1), executing a random walk as new, uncorrelated 
fluctuations are encountered. It follows that |V^/i| 2 re- 
mains smaller than e for at least a distance of order e _1 . 
(This is a conservative estimate because cusp formation, 
equivalent to discarding certain branches of a multival- 
ued eikonal solution, can and docs eliminate relatively 
large tilts.) 

But the rescaling of xm performed in Sec. IIIBI shows 
that the characteristic distance for front equilibration is 
of order e -2 / 3 <^ e -1 . Thus, at a minimum, our approxi- 
mation | V_l/i| 2 <C e remains valid well after a statistically 
steady state is reached, and its validity is then assured 
forever. Nonetheless, the contribution of |Vj_/i| 2 must 
be included in a useful reduced equation for propagation, 
because this nonlinear term is responsible for producing 
the steady state. The leading terms in Eq. (|A2[) then give 
the inviscid KPZ equation 

|L = I|V^-<x; (A3) 

the omitted terms are negligible compared to |Vj_/i| 2 . 
We conclude that the formation and properties of the 
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Huygens-propagation steady state (for sufficiently small 
e) are accurately described by the KPZ equation with 
a non- white-noise perturbation —<r and with viscosity 
taken to zero. 



so that 



A(u) = iuz(u) — ~ / dv z(v) 



(B6) 



APPENDIX B: GAUSSIAN TRIAL FUNCTIONS 
AND REPLICA SYMMETRY BREAKING 

To allow the expectation value (ip\H n \ip) of the n- 
particle Hamiltonian to be expressed analytically 
in n, we adopt the usual isotropic Gaussian trial wave 
functions 



■0(yi 



oc exp 



a, 6=1 



)ab Ya ■ Yb 



(Bl) 



which obey (^>| y a • y b \ip) = (d - l)Q a b, where Q is a 
positive-definite matrix. For positive integer n, the opti- 
mal choice of Q to approximate the ground state would 
be invariant under arbitrary permutations of the repli- 
cas, because this is a symmetry of the Hamiltonian. For 
n — > 0, however, at least within this variational approach, 
the permutation invariance is violated through hierarchi- 
cal replica symmetry breaking [28| . The hierarchical 
matrix Q is constructed as follows, starting from positive 
integer n. Given integers 1 = too < m i < • • • < n^K < 
ttik+i = n such that each rrii divides TOi+i, define Mj as 
the nxn block-diagonal matrix consisting of submatrices 
of size mi x rrii with all entries 1 (e.g., Mq is the identity 
matrix). Then define 



Q 



K+l 

h ^ 



(B2) 



i=0 



where the scalars hi have dimensions L 2 . The n particles 
are thus divided into blocks, sub-blocks, etc., that are 
bound on different length scales; this already suggests a 
connection to renormalization ideas. 

Because all the Mi are seen to commute, the eigenval- 
ues of Q are readily found. Each all-1 submatrix of Mi 
has one eigenvalue rrii and — 1 eigenvalues 0, so Mj 
has n/rrii eigenvalues rrii and n — n/rrii eigenvalues 0. For 
< i < K, the matrix Q has n/rrii — n/rrii + i eigenvec- 
tors that are in the null space of Mj for j > i but in the 
nonzero eigenspace of Mj for j < i. Thus an eigenvalue 
of Q with multiplicity n/rrii — njrrii + \ is 



(B3) 



It is convenient to define piecewise constant functions on 
the real interval 1 < u < n: 



A(u) = Ai 

i 



(rrii < u < m i+ i), (B4) 



(rrii < u < m i+ i), 



(B5) 



There is a single further eigenvalue of Q given by Eq. 
(|B3[) with i — K + l, corresponding to the eigenvec- 
tor (1, . . . , 1). This represents a center-of-mass transla- 
tion, and the eigenvalue A^+i should go to infinity in the 
ground state (complete freedom of the center of mass). 

We can interpret z(rrii) — 2j^-_ bj as the variance of 
each component of intcrparticle separation, 

z{tm) = {ip\ [n-(y a -y 6 )] 2 \iJj) = Q a a + Qb b -2Q ab , (B7) 

for indices a and b that are in the same size-mi block but 
different size-TOj_i blocks. [We see this because, from 

Eq. (JB2J, Q aa = Qbb = Ef=0 lb 3 and Qab = b 3 ■] 

Using the Gaussian identity (fl7|) . this time with ( = 
• (y a — y;,), it follows that 

exp[zk ± • (y Q - y b )] = exp[-±|kj 2 z(m t )}. (B8) 
Then, from Eq. pp|) . we obtain 
(i>\V(\y a -y b \)W 



d d ~ l \s. 



(27T) 

= B (z{mi)), 



_L e xp[-i|k ± | 2 zK)]^(|k ± |) 



(B9) 



where the function Bq is determined by the spectrum of 
the random medium. The number of index pairs (a, b) of 
the type considered is n(rrii — to^-i); the range 1 < i < 
K+l covers all pairs with a ^ b. There are an additional 
n self-pairs (a = b) for which, in place of Eq. (|B9[) . we 
have (tp\ V(0) = B (0). Thus 

n 

EM y (iy<*-y<w> 



a,b=l 



K+l 



-Bo(O) + ^2 n ( m i ~ TO i-i) B (z(rrii)) 



i=l 

" { B o (0) + duB (z(u)) 



(BIO) 



Meanwhile, for the kinetic part of the Hamiltonian, we 
find 



^iV^H-Kd-ljtrQ- 1 



A" 



n \ 1 



\{d-l)n 



m l m i+1 J A 
" du 



! u 2 k{uY 



(Bll) 
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where we have set I/Ak+i — 0. Combining these results, 
we obtain the expectation value of the Hamiltonian (|23[) . 

{mn\^) = \r\d~l)n / 

s J 1 u l A(u) 

-^{ B o(°)+J i duB (z(u))y (B12) 

This expression is very similar to the one obtained in 
a replica treatment of directed polymers focused on the 
Gaussian medium There, however, the factor i mul- 
tiplying the B terms was incorrectly omitted. Also, the 
self-interaction term cx Bq(Q) was not shown (a constant 
offset that does not affect the variational optimization 
but is important for absolute energy values); the fac- 
tor 1/r in Eq. ((24|) was included in the definition of the 
Hamiltonian H n ; and the following notations were used: 
N = d — I, /3 = 1/r, X(u) = [2A(u)}-\ Q{u) = z(u), 
f = -B /(d-l)._ 

In the formal limit n — > 0, the block "sizes" m, are 
assumed to be real numbers in the reversed sequence 
n = rriK+i < mpc < ■ ■ ■ < m\ < ttiq = 1, with arbitrar- 
ily large K and no divisibility constraints. As ^ oo, 
then, A(u) and z(u) become general functions on the in- 
terval < u < 1. Because A (it) represents an eigenvalue 
of Q and z(u) represents a variance, both must be non- 
negative. A further constraint arises when we require 
nonnegativity of variances involving arbitrary numbers 
of particles from various blocks (arbitrary because, upon 
analytic continuation in n, there is no limit on the num- 
ber of particle indices that can be formally considered, 
unlike the case of positive integer n). The parameters 
bi, which determine the matrix elements of Q, must be 
such that all eigenvalues A(it) are nonnegative, not just 
for n — > but also for arbitrary realizable integer values 
1 = mo < ■ ■ • < ttik+1 = n. From Eq. (|B3|) . we see 
that nonnegativity of A^ for arbitrarily large mi requires 
bi > 0. Consequently z(u), as defined in Eq. (|B5|) . is 
a nondecreasing function for 1 < u < n; and so in the 
ri — > limit, where the order of the rrii is reversed, z(u) 
must be a nonincreasing function for < u < 1. In fact, 
from Eq. (|B6[) and its implication 

A'(tt) = \uz'{u), (B13) 

a nonnegative and nonincreasing function z (it) will pro- 
duce a function A(w) with the same properties. 

APPENDIX C: NUMERICAL METHOD FOR 
THE INVISCID KPZ EQUATION 

We simulate the one-dimensional inviscid KPZ equa- 
tion (|85p by a Lagrangian finite-element method in which 
the only approximation, aside from the periodic domain, 
is a time and space discretization of the white noise rj. 
Equation ([85]) itself is solved exactly (except for roundoff 
error). We solve this equation, rather than just the Burg- 
ers equation obtained from it, because h remains continu- 
ous at shocks and can be used to track them, and because 



computing h allows use of the formula A = (dh/dt) (in 
which averaging over t is particularly convenient). The 
time discretization of 77 is in a sense the simplest possible: 
a sequence of "delta-function kicks" at equally spaced in- 
stants tk — kb with k = 1, 2, . . . . Each kick has a random 
x profile (of a form to be described) and an overall am- 
plitude that scales with 6 1 / 2 , to produce white noise as 
6-> 0. 

The steps in our numerical method are displayed in 
Fig. We represent a "snapshot" of h(x) by a piece- 
wise quadratic function on various x intervals (elements) . 
Continuity of h is required, but dh/dx can jump dis- 
continuously upward at shocks (corresponding to Huy- 
gens cusps that are concave, not convex). Between kicks, 
the elements evolve dynamically in a way that preserves 
the piecewise quadratic form. Specifically, the nonlin- 
ear term ^(dh/dx) 2 in Eq. ([85]) is quadratic in x if h 
is, and so when 77 = the exact solution remains in the 
piecewise quadratic space. Element boundaries without 
shocks are simply advected at the local Burgers veloc- 
ity w = —dh/dx (which is constant along Lagrangian 
trajectories) pjj. For boundaries with shocks, the adja- 
cent segments are at first allowed to overlap. Even with 
no shocks initially, elements can overlap if boundaries 
pass through one another, indicating formation of a new 
shock. 

In this way h(x) can be evolved directly to the time of 
the next kick, but it generally becomes multivalued, con- 
sisting of quadratic functions on overlapping x intervals. 
The solution is then "trimmed" to obtain nonoverlapping 
elements containing the largest value of h at each x , as 
illustrated in Fig. [2](a). Some elements are cut down; 
others are discarded entirely (such as those that turn 
inside out when a new shock forms). The portions re- 
moved can be interpreted as Burgers fluid elements that 
have run into shocks or as segments of the Huygens front 
that have been overtaken at cusps. The trimming proce- 
dure is tractable if no overlaps occur between elements 
that were more distant than next-to- nearest neighbors. 
If this condition is violated, we split the time interval 
in half and perform the evolution in two stages, recur- 
sively. Because the configuration of h(x) results from a 
random process, the next-to-nearest interaction is in fact 
sufficient for evolution over a finite time interval. That 
is, barring exact synchronization between different parts 
of space (which is vanishingly unlikely), any overlap of 
more distant elements can be reduced to discrete stages 
of evolution in which only nearest and next-to-nearest 
neighbors overlap. For example, a shock gradually ab- 
sorbs the elements on either side of it (say elements 1 
and 2), but one of the two will disappear first, result- 
ing in a renumbering of the remaining elements. Only if 
they disappeared at the same instant would we have an 
unavoidable interaction between elements and 3. 

"Kicking," illustrated in Fig. [2jb), preserves the as- 
sumed form of h(x) if the spatial profile of the kick is 
also piecewise quadratic. Such a profile is synthesized 
from a given spectrum D(k) by first generating delta- 
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FIG. 2: Lagrangian finite-element method for the inviscid KPZ equation. Plots show stages in the simulated evolution of h(x) 
on a periodic x domain, (a) Filled circles bound quadratic elements of the initial configuration, distinguished by alternating 
dotted, dashed, and solid curves. One initial cusp is present, in the right half of the plot. After advection of all boundaries, 
h(x) consists of corresponding elements bounded by open circles, and is no longer single- valued. The filled circle at the initial 
cusp separates into two open circles, and a dashed element in the left half of the plot turns inside out when its boundaries pass 
through one another (inset). The trimming procedure then constructs a single- valued h{x) by retaining only the largest values 
(shaded band). As a result, a solid element disappears into the existing cusp, and the inside-out dashed element disappears to 
form a new cusp, (b) Open circles indicate trimmed element boundaries from (a), including updated cusp locations. Vertical 
lines form grid for kicking. Arrows indicate deformation of kicking boundaries to existing element boundaries, except third grid 
line from left, which has no nearby element boundary and introduces a new one. Kicking produces elements bounded by filled 
circles, without altering the number or location of cusps. This final h(x) can then be evolved again as in (a). 



function spikes on a uniform x grid (spacing 5) based on 
a spectrum D{k), and then forming a smooth piece- 
wise quadratic function by in effect analytically integrat- 
ing three times with respect to x. (This kick profile has 
an everywhere continuous derivative to avoid introduc- 
ing additional shocks, especially ones of unphysical sign.) 
The choice of the forcing grid is the only way in which 
a finite spatial resolution 6 enters the simulation. As 
described, the kicking process would introduce a new el- 
ement boundary in h(x) at every grid point, since gener- 
ically every boundary introduced in a previous kick has 
moved at least slightly. Each boundary survives for some 
characteristic timescale before vanishing into a shock, 
and so the steady-state average number of boundaries 
present is proportional to the rate at which they are intro- 
duced, which would scale with the kicking rate 1/6. As we 
take b — > to represent white noise accurately, we would 
have an explosion in the number of elements, but they 
would be mostly redundant since the spatial resolution 
of the driving noise is still limited by the grid. We adopt 
a more efficient approach that adds a new boundary only 
when h{x) is insufficiently resolved in the neighborhood 
of the point in question, and otherwise deforms the kick 
so as to take advantage of an existing nearby boundary 
rather than insisting on the planned grid. The small- 
scale deformation distorts the high- wave-number part of 
the noise spectrum and thus requires a somewhat finer 
grid to achieve the same accuracy. But the cost is out- 



weighed by the much-improved behavior as b — > 0: Now 
the average number of boundaries remains fixed. 

Conventional numerical methods for the Burgers-KPZ 
and similar equations require fine spatial grids to resolve 
shocks, and nonzero viscosity to stabilize them. By con- 
trast, we take the inviscid limit from the start, allowing 
an explicit geometric representation of shocks. The spa- 
tial grid need only resolve the forcing; inviscid shocks 
are perfectly sharp and do not introduce smaller length 
scales. If the forcing is spatially smooth (as for media 
G and MG), then there is a significant advantage in ef- 
ficiency from using a relatively coarse grid and still cap- 
turing sharp shocks. With spatially rough forcing (i.e., a 
spectrum with a long high-wave-number tail, as for media 
E and ME), the advantage is less clear because a fine grid 
is needed anyway to resolve the forcing accurately. Nev- 
ertheless, because we work directly at Rcb = oo, there 
is one fewer parameter contributing systematic errors. 
Our method is designed explicitly for a one-dimensional 
simulation (d = 2); generalization to the Burgers-KPZ 
equation in two or more dimensions (d > 3) is possible 
in principle, but the computational geometry would be 
much more intricate and possibly intractable. As with 
conventional methods, the computational cost of simula- 
tions in higher dimensions would be severe. 

We now describe how the parameters of our simula- 
tions are chosen and how the uncertainties in the results 
are estimated. By systematic error we mean the differ- 
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ence between the precise average of a quantity over many 
runs of a practical simulation (where each run takes a fi- 
nite computation time) and the precise average for the 
idealized problem statement (where space and time are 
considered infinite and continuous). An efficient compu- 
tational approach involves balancing statistical and sys- 
tematic errors, both of which contribute to the overall 
uncertainty of the result. In our case, the relevant sys- 
tematic parameters that ideally approach infinity are L, 
1/(5, 1/6, and the time T allowed for equilibration before 
data are taken. We estimate A by subsequently averag- 
ing dh/dt over a time interval 3T, based on the follow- 
ing considerations: It cannot be efficient to spend much 
more time on equilibration than on averaging, and even 
if it were ideal to spend only a tiny fraction of the com- 
putation time on equilibration, spending j on it reduces 
the available data for averaging only modestly, increasing 
the statistical error by (I) 1 / 2 — 1 — 15%. 

Our framework for treating systematic errors is a con- 
servative assumption that these errors scale with the re- 
ciprocal of the parameters given above. (This is the slow- 
est convergence that we would reasonably anticipate.) 
That is, we assume that the precise average computed 
from many runs of a simulation tends to the true value 
of A with asymptotic corrections of order l/L, order 6, 
order &, and order 1/T. Consequently, we can extrapo- 
late from simulations performed with different finite pa- 
rameters to estimate the result of an ideal simulation. 
Because we do not trust this extrapolation as a quanti- 
tative model, we will apply it only when all the simula- 
tion results contributing to the extrapolation are statis- 
tically indistinguishable (consistent with identical under- 
lying averages). 

Specifically, our final extrapolation will be based on 
a "most refined" simulation, with parameter set a, and 
N = 4 lesser simulations /?o, • ■ ■ , Pn-i, each based on 
a with one parameter halved. All these simulations are 
repeated as necessary to obtain averages A and Bi with 
some statistical uncertainty a (to be determined below) . 
Because of the parameter halving and the assumed recip- 
rocal scaling of systematic errors, the extrapolation takes 
the simple form 

A~ A+(A-B ) + --- + (A-B N - 1 ). (CI) 

To ensure that the result is fairly insensitive to our spe- 
cific model of systematic errors, we require each correc- 



tion term A — Bi to be within two standard deviations 
of zero, i.e., 

\A- Bi\ < 2V2a. (C2) 

(Taking a difference of two independent random variables 
multiplies the standard deviation by y/2.) If the decay of 
systematic errors is more rapid than assumed, Eq. (|C1[) 
is needlessly imprecise but not incorrect, because the ex- 
aggerated correction terms are statistically equivalent to 
zero. Our final estimate of A, being a sum of indepen- 
dent random variables (N + 1)A — Bq — ■ ■ ■ — Bjv-i, has 



TABLE III: Most refined numerical simulation parameters, 
and resulting extrapolated estimate of A, for media in Table 
E(with a = 1). 



Medium 


log 2 i 


log 2 (3<5) 


log 2 b 


log 2 T 


A 


G 


8 


-2 


-4 


8 


1.535(14) 


MG 


7 


-3 


-6 


6 


1.450(13) 


E 


7 


-5 


-5 


5 


1.73(4) 


ME 


6 


-5 


-5 


4 


1.66(3) 


a variance 


[(N + l) 2 + N]a 2 


= 29a 2 . 


Thus a 


should be 



taken as the desired overall uncertainty divided by a/29- 
By repeated doubling, we can locate a parameter set a 
sufficiently refined that Eq. (|C2[) holds. 

The purpose of the extrapolation method is to obtain 
a conservative assessment of the uncertainty contributed 
by systematic errors. Our procedure is qualitatively sim- 
ilar to traditional rules of thumb, such as, "The change in 
the simulation result upon doubling the resolution is an 
estimate of the systematic uncertainty." But the some- 
what arbitrary technique of doubling the resolution (or 
other parameters) is replaced here by a more fundamen- 
tal assumption about the scaling of systematic errors. (In 
simpler problems, such as the finite-difference solution of 
deterministic differential equations, the correct scaling 
can be readily obtained by analysis, but we desire a ro- 
bust "black box" method.) 

For each of our four example spectra, Table IIIII gives 
the parameters a of the most refined simulation per- 
formed, and the resulting extrapolated estimate of A in- 
cluding both statistical and systematic uncertainties. 
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